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EVALUATION  OF  ATTENUATION/MINIMUM-PHASE  PAIRS 
BY  MEANS  OF  TWO  FAST  FOURIER  TRANSFORMS 

INTRODUCTION 

It  is  often  important  to  determine  whether  a  given  linear 
device  is  minimum-phase  [1],  because  if  so,  it  is  then  possible 
to  compensate  the  filter  characteristic  with  reciprocal  pole-zero 
locations  and  obtain  an  overall  all-pass  characteristic  with  flat 
amplitude  and  linear  phase  responses.  A  relatively  simple  way  of 
making  this  determination  is  to  measure  the  attenuation  (or 
decibel  gain)  and  actual  phase  shift  of  the  given  linear  device 
and  then  compute  the  minimum-phase  corresponding  to  the  measured 
attenuation.  If  this  latter  calculated  phase  agrees  with  the 
actual  measured  phase,  then  the  filter  is  minimum-phase. 

The  minimum-phase  corresponding  to  a  given  attenuation 
function  is  determined  analytically  by  a  Hilbert  transform 
[2;  chapter  6,  article  22]  or  [3;  section  10-3].  However,  this 
direct  integral  evaluation  is  computationally  unattractive  due  to 
two  poles  on  the  line  of  integration  [3;  (10-67)].  In  addition, 
it  yields  only  a  single  value  for  the  phase  after  each  numerical 
integration.  We  will  circumvent  both  of  these  difficulties  by 
first  subtracting  the  singularities  (which  will  be  handled 
analytically)  and  then  employing  fast  Fourier  transforms  for 
efficient  numerical  evaluation  of  the  entire  phase  response. 
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TRANSFER  FUNCTION  RELATIONS 

FILTER  CHARACTERIZATIONS 

A  linear  time-invariant  filter  is  characterized  by  its 
impulse  response  h(x)  or  by  its  transfer  function  H(f)  according 
to  Fourier  transform 

H ( f )  -  J  dx  exp(-i2nfx)  h(x)  ■  F{h(x)}  .  (1) 

(Integrals  without  limits  are  over  the  range  of  nonzero 
integrand.)  Both  the  impulse  response  h(x)  and  the  transfer 
function  H(f)  can  be  complex  functions  of  time  delay  x  and 
frequency  f,  respectively. 

The  transfer  function  will  be  represented  in  terms  of  its 
real  and  imaginary  parts  according  to 

H ( f )  -  Hr(f)  +  i  H.(f)  ,  (2) 

where 

Hr(f)  -  j(H(f)  +  H* ( f )  ]  , 

Hi ( f )  -  jjlH(f)  -  H* ( f ) ]  .  (3) 

It  can  also  be  represented  in  terms  of  its  even  and  odd  parts  as 

H ( f )  -  He(f)  +  H0(f)  ,  (4) 

which  are  generally  defined  according  to 
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He(f)  “  |[H(f)  +  *  J  dT  cos(  2nfx)  h(T) 


HQ(f)  -  |[H(f)  -  H ( -f ) ]  -  -ij  dx  sin( 2 nfx)  h(x) 


(5) 


Functions  Hg(f)  and  HQ(f)  are  both  complex  generally,  whereas 
Hr(f)  and  H^(f)  are  always  real.  Impulse  response  h(x)  can  be 
complex . 

(In  the  special  case  where  impulse  response  h(x)  is  real, 

then 


HQ(f)  »  i  H^(f)  -  -ij  dx  sin(2nfx)  h(x)  .) 


(6) 


CAUSAL  FILTER 

A  filter  is  said  to  be  causal  when  its  impulse  response  h(x) 
is  zero  for  negative  arguments;  that  is. 


h(x)  -  0  for  x  <  0  . 


(7) 


However,  h(x)  can  still  be  a  complex  function  of  x.  In  this 
causal  case,  the  real  and  imaginary  parts  of  the  transfer 
function  H(f)  satisfy  a  pair  of  Hilbert  transform  relationships, 
provided  that  h(x)  does  not  contain  any  impulses  at  the  origin; 
see  also  [3;  page  198].  The  Hilbert  transform  of  an  arbitrary 
complex  function  G(x)  is  defined  as 
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s'0'*"  *  H du  ^  ■  «  ®  Gu>  •  <8) 

where  the  tic  mark  on  the  integral  sign  denotes  a  principal  value 
integral  [4;  section  3.05]  and  9  denotes  convolution.  Principal 
value  integrals  are  considered  in  appendix  A. 

In  order  to  derive  the  Hilbert  relations  of  interest,  let 
U(x)  be  the  unit  step  function, 

f 1  for  x  >  O') 

U(x)  -  {  [  .  (9) 

lo  for  x  <  oj 

Then,  because  h(x)  is  causal,  transfer  function  (1)  becomes 

H ( f )  -  J  dx  exp(-i2nfx)  h(x)  U(x)  »  P{h(x)  U(x)}  - 

-  F{h( x) }  9  F{U(  x)  }  -  H  (  f )  9  [|i(f)  +  - 

-  |  H ( f )  -  |  H  { H  (  f  )  }  .  (10) 

Here,  we  used  the  Fourier  transform  of  the  unit  step  function 
U(x)  [3;  (3-13)]  and  definition  (8).  Equation  (10)  yields 

H(f)  -  -i  H{ H( f ) }  (11) 

or,  more  explicitly, 

Hr(f )  -  H{H.(f ) }  -  •  Hi(f)  , 

H.(f)  -  -  H{Hr(f)}  -  -  ^  ©  H  r  (  f )  .  (12) 
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We  repeat  that  transfer  function  relations  (12)  hold  true  even 
when  impulse  response  h(T)  is  complex;  only  causality  is  used. 
Analogous  properties  to  (12)  hold  between  the  even  and  odd  parts, 
Hg(f)  and  H  (f),  of  the  transfer  function  H(f)  as  well.  Namely, 
because  the  Hilbert  transform  of  an  even  (odd)  function  is  odd 
(even),  there  follows,  for  a  causal  (but  possibly  complex)  h(x), 

He(f)  -  -i  H{HQ(f)}  ,  HQ(f)  -  -i  H{He(f)}  .  (13) 

If  h(x)  contains  an  impulse  at  the  origin,  both  parts  of 
(12)  are  false,  even  though  h(x)  may  be  causal.  Consider 

h(x)  -  (a  +  ib)  6(x),  a  and  b  real  .  (14) 

Then  (1)  yields  constant  transfer  function 

H ( f )  -  a+ib,  Hr(f)  -  a,  H± ( f )  -  b,  H0(f)  -  a+ib,  HQ(f)  -  0.  (15) 

But  since  the  Hilbert  transform  of  a  constant  is  zero 

[4;  section  3.05),  neither  part  of  (12)  is  satisfied,  and  the 

first  part  of  (13)  is  false.  On  the  other  hand,  if 

h(x)  -  (a  +  ib)  5(x  -  T)  ,  a  and  b  real  ,  (16) 

then  (12)  and  (13)  are  satisfied  only  if  T  >  0.  Here,  we  used 
the  facts  that 

H{cos(2nfT))  -  sin(2nf|T|),  H{ sin( 2 nfT) }  -  -sgn(T)  cos( 2 nfT ) , ( 17 ) 

where  sgn(T)  is  the  polarity  of  T.  Henceforth,  we  assume  that 
components  like  (14)  and  (15)  are  not  present  in  the  filters  of 
interest;  see  also  [3;  page  198). 
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For  a  causal  filter,  (2)  and  (12)  afford  a  method  of 
obtaining  the  complete  transfer  function  from  its  real  part 
alone,  according  to 


H ( f )  -  Hr(f)  +  i  H.(f)  - 

-  Hr(f)  -  i  H{Hr(f)}  .  (18) 

However,  a  more  attractive  approach,  computationally,  is  to  use 
Fourier  transforms,  as  follows.  Define  inverse  Fourier  transform 

h(r)  «  F-1 { H r ( f ) )  -  J  df  exp( i2nf t)  Hr ( f )  (19) 

for  any  real  part  Hr(f).  (The  notation  hr(t)  cannot  be  used 
instead  of  h(x),  because  h(t)  is  not  the  real  part  of  h(t),  nor 
is  h(x)  necessarily  real.)  Substitution  of  (3)  into  (19) 
immediately  yields 

h(x)  -  |[h(T)  +  h*(-x)]  ;  h(-x)  -  h*(x)  .  (20) 

(These  particular  relations  in  (20)  actually  hold  true  for  any 
filter  h(t),  noncausal  as  well  as  complex.)  Then  because  h(x)  is 
causal,  there  follows  directly 

f2h(x)  for  t  >  O') 

h(T)  -  {  \  -  2  h(x)  U(t)  .  (21) 

l  0  for  t  <  Oj 

In  summary,  the  method  for  obtaining  the  complete  transfer 
function  H(f)  from  just  its  real  part  Hf(f),  for  a  causal  filter, 
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is  to  perforin,  in  order,  the  following  operations: 

h(x)  -  F_1 { Hr ( f ) }  , 
h(x)  -  2  h  ( x )  U(x)  , 

H ( f )  «  h ( x ) }  .  (22) 

This  procedure  requires  two  Fourier  transforms,  which  can  be 
accomplished  very  quickly  and  efficiently  by  means  of  two  fast 
Fourier  transforms.  Furthermore,  a  fast  Fourier  transform 
output  sweeps  out  the  complete  range  of  argument  values,  whereas 
the  brute  force  Hilbert  transform  integral  of  (18)  and  (8) 
requires  an  additional  numerical  integration  for  each  frequency  f 
of  interest.  Functions  h(x)  and  h(x)  in  (22)  can  be  complex. 

An  accuracy  check  on  the  procedure  in  (22)  is  afforded  by 
comparing  the  real  part  output  of  the  Fourier  transform  in  the 
bottom  line  with  the  input  Hf(f)  utilized  in  the  top  line.  The 
complete  set  of  function  values  of  H^ff)  for  all  f  is  required 
for  this  procedure;  in  return,  the  complete  set  of  values  of 
H^(f),  for  all  f,  results.  The  operations  in  (22)  are  linear 
insofar  as  the  overall  transformation  of  Hf(f)  is  concerned,  and 
so  superposition  can  be  used  for  any  breakdown  of  Hf(f)  into 
components,  if  desired. 

The  rule  for  obtaining  H(f)  or  H^(f)  from  H  (f),  as  given  in 
(22),  applies  whether  filter  H(f)  is  minimum-phase  [1]  or  not. 

The  only  prerequisite  for  the  validity  of  (22)  is  the  causality 
of  impulse  response  h(x). 

If  only  H  (f)  were  available  (instead  of  H  ( f ) ) ,  a  more 
e  r 

attractive  procedure  for  obtaining  H(f)  or  H^(f)  than  using  (4) 
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and  Hilbert  transform  (13),  is  to  observe  that,  in  general,  for 
any  filter,  the  inverse  Fourier  transform 

F_1{He(f)}  -  J  df  exp( i2 nf t)  He(f)  -  ^[h(x)+h(-x) ]  -  hft(x).  (23) 

Here,  we  used  (5),  the  inverse  to  (1),  and  the  general  definition 
of  the  even  part  of  an  arbitrary  complex  function.  Then,  if  h(x) 
is  causal,  we  have 

h(x)  -  2  he(x)  U ( t )  .  (24) 

Thus,  the  procedure  for  obtaining  H(f)  is  identical  to  (22)  if  we 
replace  H r ( f )  and  h(x)  by  H& ( f )  and  hg(x),  respectively. 

ONE-SIDED  SPECTRAL  FUNCTIONS 

The  analogous  situation  in  the  frequency  domain  (to  causality 
in  the  time  delay  domain)  is  as  follows:  if  (complex  function) 
A(f)  is  zero  for  negative  arguments,  that  is, 

A(f)  -  0  for  f  <  0  ,  (25) 

then  a  procedure  similar  to  (10)— (11)  reveals  that  the  inverse 
Fourier  transform  of  A(f)  is  given  by 

a(x)  ■  r"1{A(f)}  -  i  H{a( x) }  .  (26) 

That  is,  in  terms  of  real  and  imaginary  parts, 

af(x)  -  -  H{ai(x)}  ,  ai(x)  -  H{ar(x))  .  (27) 

The  function  a(x)  is  called  an  analytic  waveform,  for  reasons  to 
become  apparent  shortly. 
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GENERAL  SPECTRAL  RELATIONS 

For  future  purposes,  the  Hilbert  transform  of  a  completely 
arbitrary  complex  waveform  b(x), 

bH(r)  *  H { b ( t ) }  -  ©  b(T)  ,  (28) 

n  —  JIT 

has  spectrum  (Fourier  transform) 

f-i  B(f)  for  f  >  O') 

F{b„(x)}  -  -i  sgn(f)  B(f)  -  {  l  ,  (29) 

l  i  B(f)  for  f  <  Oj 

where  B(f)  is  the  spectrum  of  b(x).  Here,  we  used  the  fact  that 
the  following  two  functions  are  a  Fourier  transform  pair 
(3;  apply  (2-34)  to  (3-9)]: 

~  -i  sgn(f )  .  (30) 

The  left-hand  side  of  (29)  is  the  Fourier  transform  of  the 

Hilbert  transform  of  b(x).  It  cannot  be  labeled  as  B„(f),  which 

n 

is  the  Hilbert  transform  of  the  Fourier  transform  B(f)  of  b(x). 
The  two  operations  of  Hilbert  transformation  and  Fourier 
transformation  are  not  interchangeable,  in  general. 

It  follows  from  (29)  that 

F{ b ( x )  +  i  bH(x)}  -  2  B ( f )  U(f)  ,  (31) 

which  is  a  one-sided  spectrum.  Also,  b(x)  +  i  bH(x)  is  an 
analytic  waveform.  Waveform  b(x)  is  completely  arbitrary  here. 
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ANALYTICITY  OF  TRANSFER  FUNCTION 


Consider  the  causal  exponential  impulse  response 


h(x)  ■  exp(-x)  U(t)  . 


(32) 


The  corresponding  transfer  function  is 


which  has  a  pole  in  the  upper-half  f-plane  at  f  -  i/(2n),  but 
which  is  analytic  in  the  lower-half  f-plane.  (The  lower-half 
f-plane  corresponds  to  the  right-half  s-plane  of  Laplace 
transforms . ) 

This  analyticity  of  the  transfer  function  H(f)  in  the  lower- 
half  f-plane  is  generally  true  for  causal  finite-energy  filters, 
as  may  be  seen  by  the  following  argument.  Let  frequency  f  be  a 
complex  variable  with  real  and  imaginary  parts  according  to 
f  -  f  +  if..  Then,  for  a  causal  filter,  (1)  can  be  expressed 
more  explicitly  as 


H  ( f ) 


(34) 


0 


The  first  exponential  in  (34)  has  magnitude  1  for  all  t  on  the 
contour  of  integration.  And  if  f^  <  0,  the  second  exponential 
term  in  (34)  decays  with  increasing  r,  keeping  the  integral 
convergent,  as  it  was  for  f^  -  0.  That  is,  transfer  function 
H(f)  is  analytic  in  the  lower-half  f-plane  for  a  causal  impulse 
response  h(x).  Notice,  however,  that  no  statements  can  be  made 
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about  the  locations  of  the  zeros  of  transfer  function  H(f)  in  the 
complex  f-plane.  Thus  we  have 

causal  h(t)  — >  analytic  H(f)  in  lower-half  f-plane  .  (35) 

The  converse  is  also  true,  namely,  that  analyticity  implies 
causality.  To  develop  this  point,  express  the  inverse  Fourier 
transform  to  (1)  in  the  form 


i ( t )  -  J  df  exp( i2nf t)  H(f) 


J  df  exp(i2nfrT)  exp(-2nf^r)  H(f)  ,  (36) 


where  contours  and  C2  are  depicted  in  the  complex  f-plane  in 
figure  1.  Because  transfer  function  H(f)  is  analytic  in  the 
( crosshatched )  region  between  contours  and  C2 ,  we  are  allowed 
to  move  the  integration  freely  between  them,  as  done  in  (36), 


Figure  1.  Complex  f-Plane  Contours 
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without  altering  the  value  h(x)  of  the  integral.  On  contour  C2* 
we  have  f^  <  0  everywhere.  Therefore,  if  t  <  0  in  (36),  the 
second  exponential  decays  to  zero  as  contour  C2  is  moved  farther 
down  in  the  f-plane.  Because  H(f)  is  analytic  in  the  lower-half 
f-plane,  we  can  move  C2  arbitrarily  far  down,  causing  the 
integrand  of  (36)  to  go  to  zero,  thereby  leading  to  a  zero  value 
for  h(x)  whenever  x  <  0.  Thus,  we  have 


analytic  H(f)  in  lower-half  f-plane  — causal  h(x)  .  (37) 


This  equation  is  the  converse  to  (35). 

Because  we  have  already  shown  in  (10)— (12)  that  a  causal 
impulse  response  h(x)  leads  to  a  transfer  function  H(f)  with 
Hilbert  transform  relations  between  its  real  and  imaginary  parts, 
it  follows  from  (37)  that  an  analytic  transfer  function  H(f) 
leads  to  the  same  conclusions.  This  means  that,  for  an  analytic 
transfer  function  H(f)  in  the  lower-half  f-plane,  we  can  use  the 
efficient  procedure  given  in  (22),  in  terms  of  two  (fast)  Fourier 
transforms,  to  find  the  imaginary  part  H^(f),  given  only  the  real 
part  Hr( f ) . 

For  the  example  given  earlier  in  (33),  we  have  real  part 


Hr(f ) 


1 

1  +  ( 2nf ) 2 


Then  from  (22),  we  obtain,  in  order, 


h(x)  -  j  exp(-|x|)  ,  h(x)  -  exp(-x)  U(x)  ,  H(f)  -  ^  -  *2nf  * 
which  corroborates  (32)  and  (33). 
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MINIMUM-PHASE  TRANSFER  FUNCTIONS 

From  this  point  on,  we  presume  that  impulse  response  h(x)  is 
causal  and  that  transfer  function  H(f)  contains  only  poles  and 
zeros.  It  then  follows  from  (35)  that  transfer  function  H(f) 
has  no  poles  in  the  lower-half  f-plane.  We  also  assume  now  that 
H(f)  has  no  zeros  in  the  lower-half  f-plane;  that  is,  the  filter 
is  minimum-phase  [1,2,3].  In  this  case,  the  function 

Q(f)  -  -  In  H(f)  (38) 

is  analytic  in  the  lower-half  f-plane,  because  the  function  In  z 
is  nonanalytic  only  at  z  «  0  and  z  «  •  in  the  complex  z-plane. 


Accordingly,  by  analogy  to  (37),  inverse  Fourier  transform 

q(x)  -  J  df  exp(i2nfx)  Q(f)  (39) 

is  causal.  (An  example  is  given  in  appendix  B.)  Therefore,  just 
as  shown  in  (10)  —  (12) ,  the  real  and  imaginary  parts  of  Q(f), 

Q(f)  -  Qr(f)  +  i  Qj ( f )  ,  (40) 

can  be  found  from  each  other  by  means  of  Hilbert  transforms.  In 
particular,  as  in  (12), 

Qr(f)  -  H{Q.(f)}  ,  Q.(f)  -  -  H{Qr(f)}  .  (41) 


Alternatively,  according  to  the  sequel  to  (37),  because  Q(f) 
is  analytic  in  the  lower-half  f-plane,  the  imaginary  part  Q^(f) 
can  be  found  from  real  part  Qr(f)  according  to  procedure  (22) 
involving  two  Fourier  transforms. 
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Interesting  interpretations  of  minimum-phase  filters,  in 
terms  of  their  group  delay  and  rate  of  energy  flow  through  the 
filter,  are  given  in  [5;  pages  132-3].  In  particular,  the 
minimum-phase  filter  has  the  smallest  group  delay  of  any  stable 
filter  with  specified  magnitude  transfer  function. 


ATTENUATION  AND  PHASE 


There  is  another  way  of  describing  a  transfer  function  H(f) 
rather  than  by  its  real  and  imaginary  parts,  which  is  very  useful 
in  some  applications.  Namely,  let 


where 


H ( f )  -  exp[-a(f)  -  i  0(f)]  , 


a(f ) 
0(f) 


attenuation 
phase  shift 


of 


filter  . 


(42) 


(43) 


Reference  to  (38)  and  (40)  immediately  reveals  that 

a(f)  -  Qr(f)  ,  0(f)  -  Qi(f)  .  (44) 

Therefore,  if  filter  H(f)  is  minimum-phase,  according  to  the 
discussion  in  (38)-(41),  a(f)  and  0(f)  can  be  found  from  each 
other  by  means  of  Hilbert  transforms.  In  particular, 

0(f)  -  -  H{  a(  f )  }  -  -  -^1  ©  a(f)  .  (45) 

(Strictly,  this  relation  is  not  usable  and  must  be  modified  to 
allow  for  attenuations  a(f)  with  logarithmic  singularities;  for 
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example,  see  [3;  pages  206-8].  This  manipulation  is  discussed 
in  appendix  C . ) 

Alternatively,  the  procedure  in  (22)  can  be  employed  in  the 

form 

3(t)  *  F— 1 { o( f ) }  , 
q(x)  -  2  3(t)  U(t)  , 

«(f)  +  i  3(f)  -  F{q(x)}  .  (46) 

The  function  <j(x)  is  defined  by  the  inverse  Fourier  transform  in 
the  top  line  of  (46).  Phase  shift  8(f)  for  a  minimum-phase 
filter  is  given  by  the  imaginary  part  of  the  Fourier  transform  in 
the  bottom  line  of  (46). 

A  common  alternative  descriptor  of  the  frequency  behavior  of 
a  filter  is  the  gain  G(f)  in  decibels,  defined  as 

G(f )  -  20  log10  | H ( f ) |  .  (47) 

Because  the  attenuation  follows  from  (42)  as 


ct(f)  -  -  In  t  H  (  f  )  |  , 


(48) 


the  gain  G(f)  and  the  attenuation  a(f)  are  related  by 


G(f ) 


20 


ln(10) 


o(f)  -  -  8.686  a(f) 


(49) 


Measurement  of  either  one  is  sufficient  to  find  the  other  and  to 
thereby  determine  the  phase  shift  8(f)  of  a  minimum-phase  filter. 
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EXAMPLE  AND  LIMITATION 

We  again  consider  the  example  given  in  (32)— (33),  namely 

h(  t)  -  exp(-x)  U(  x )  ,  H  ( f )  -  Y~~^i2nf  *  <5°) 

The  attenuation  and  phase  follow  from  (42)  according  to 

ct(f)  -  |  ln(  1  +  4n2f2 )  , 

(3(f)  -  arctan(2nf)  .  (51) 

If  we  attempt  to  apply  the  inverse  Fourier  transform  in  the  top 
line  of  (46)  to  the  attenuation  a(f)  in  (51),  we  encounter  a 
divergent  integral  because  o(f)  ~  ln|f|  as  f  -»  +«. 

More  generally,  if  filter  H(f)  has  a  zero  at  a  frequency  f 
equal  to  any  finite  real  value,  the  attenuation  a(f)  has  a 
logarithmic  singularity  at  that  real  frequency,  and  the  inverse 
Fourier  transform  in  (46)  diverges.  Because  typical  filters  very 
often  have  this  feature  (and  almost  always  at  f  -  0  and  f  -  +*) , 
a  way  must  be  found  to  circumvent  the  divergent  part  of  the 
inverse  Fourier  transform  integral,  so  that  the  efficient 
procedure  of  (46)  can  be  salvaged. 
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SUBTRACTION  OF  SINGULARITY 

The  procedure  to  be  used  here  is  one  commonly  adopted  to 
numerically  evaluate  convergent  integrals  with  singular 
integrands;  it  is  illustrated  by  the  example 

a 

I  -  [  dx  £212!  ,  v  <  1  .  (52) 

0  x 

If  v  is  positive,  the  integrand  has  an  infinite  cusp  at  the 
origin,  yet  the  integral  converges,  because  v  <  1.  We  express 

a  a  a 

i  -  r  dx co5x  - 1  * 1 .  r  dx cosx  — 1 1  r  dx  ^  ,  (53> 

0  x  0  x  0  xV 

which  is  allowed,  because  both  integrals  converge.  The  last 

integral  in  (53)  can  be  done  in  closed  form,  yielding  a*_v/(l-v). 
Also,  the  middle  integrand  now  behaves  as  x^-v  as  x  -»  0+,  which 
is  zero  at  the  origin,  because  2-v  >  1;  this  behavior  enables  a 
straightforward  numerical  evaluation  of  the  middle  integral. 

The  key  to  this  procedure  is  to  find  a  component  that  can  be 
integrated  in  closed  form  and  that,  when  subtracted  from  the 
given  integrand,  yields  a  well-behaved  residual  for  numerical 
integration . 


19 


TR  8667 


APPLICATION  TO  FILTERS 

The  way  we  apply  this  subtraction  procedure  to  a  given 
attenuation  a(f)  with  logarithmic  singularities  is  to  break  it 
into  two  parts, 

a(f)  «  (  f  )  +  <*2  (  f  )  ,  (54) 

where  attenuation  a^(f)  contains  all  the  singular  components  and 
has  a  known  closed  form  minimum-phase  pair  ( f ) .  (An  example  is 
furnished  by  (50)  and  (51);  some  additional  examples  are  listed 
in  appendix  D.)  Then  residual  attenuation  a2(f)  is  f°und 
according  to 

a2(f)  -  a(f)  -  a^f)  (55) 

and  is  well-behaved  for  all  f.  Residual  a2(f)  is  subjected  to 
the  repeated  Fourier  transform  procedure  detailed  in  (46), 
resulting  in  phase  shift  function  82(f)*  Finally,  the  complete 
minimum-phase  corresponding  to  the  given  attenuation  a( f )  is 
obtained  from 

8(f)  -  p1(f)  +  02(f)  .  (56) 

The  procedure  can  be  summarized  as  follows: 

a(f)  — y  8(f)  desired  ; 

a^f)  +  a2(f)  — ►  8x(f)  +  82(f)  used  .  (57) 

The  exact  choice  of  attenuation/minimum-phase  pair  a^(f), 
8^(f)  is  not  critical,  except  that  residual  <*2 ( f )  must  not  have 
any  singularities  and  must  decay  (rapidly)  to  zero  for  large  f. 
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Of  course,  the  given  attenuation  a(f)  must  be  known  for  all  f  in 
order  to  apply  this  (or  any)  procedure  for  obtaining  minimum- 
phase  shift  8(f),  whether  obtained  directly  by  Hilbert  transforms 
or  by  means  of  a  Fourier  procedure.  The  actual  numerical 
evaluation  of  the  Fourier  procedure  delineated  in  (46)  is 
accomplished  by  means  of  fast  Fourier  transforms;  the  details  are 
presented  in  appendix  E. 

SHORTCOMING  OF  HILBERT  TRANSFORM 

Suppose  that  two  minimum-phase  filters  H0 ( f )  and  H^(f)  differ 
only  by  a  complex  scale  factor: 

Hb(f)  “  c  Ha(f)  ‘  (58) 

Then 

°b(f)  "  aa(f)  '  ln|c|  ' 

8b(f)  «  8a(f)  “  arg(c)  +  2nn  ,  n  integer  .  (59) 

However,  if  a  ( f )  and  8. ( f )  are  a  Hilbert  transform  pair,  a.  (f) 
a  a  d 

and  8b(f)  cannot  possibly  be  (unless  c  *  1  and  n  -  0)  because  the 
Hilbert  transform  of  a  constant  is  zero.  Functions  a.  (f)  and 
8b(f)  are  both  "incomplete,"  in  that  attenuation  a^Cf)  contains 
no  information  about  arg(c),  while  phase  8b(f)  contains  no 
information  about  |c|.  This  means  that  the  Hilbert  transform  of 
a  given  attenuation  (phase)  yields  a  phase  (attenuation)  function 
that  can  differ  from  the  actual  phase  (attenuation)  of  a  minimum- 
phase  filter  by  an  arbitrary  additive  constant.  Some  information 
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is  inherently  absent  from  a  given  attenuation  (phase)  function. 

In  addition,  because  the  Hilbert  transform  of  a  constant  is  zero, 
additive  constants  are  lost  through  this  transformation.  (The 
situation  is  somewhat  similar  for  the  Fourier  transform  procedure 
given  in  ( 46  )  .  ) 

Alternatively,  suppose  that 


hfa(x)  -  hg(T  ~  T)  *  Hfa(f)  “  Hg(f)  exp( -i 2 Jif T )  .  (60) 


Then  filter  H^(f)  contains  a  transfer  function 
exp( -i 2 nf T ) ,  with  corresponding  attenuation  0 
Thus,  the  attenuation  contains  no  information 
delay.  However,  it  should  be  noted  that  this 
exp(-i2nfT)  does  not  possess  poles  and  zeros  a 
has  an  essential  singularity  at  f  ■  ®. 


component  of 
and  phase  2nfT. 
about  a  pure  time 
component 

t  all,  but  in  fact 
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APPLICATION  TO  MEASURED  DATA 

In  this  section,  we  will  apply  the  previous  Fourier  procedure 
to  a  measured  pair  of  attenuation  and  phase  shift  functions  in 
an  effort  to  determine  if  the  filter  is  minimum-phase.  The 
particular  filter  is  a  J15-1  transducer  used  as  a  continuous-wave 
source  in  the  10  to  900  Hertz  range.  The  transmitting  current 
response  of  this  device  is  defined  as  the  ratio 

output  pressure  (61) 

input  current 

and  is  the  transfer  function  of  interest.  The  reference  level  is 
taken  as  1  /uPa/Amp.  The  measurements  procedure  include  a 
water-path  propagation  delay  (of  unknown  value)  between  the 
transducer  and  a  calibrated  receiving  hydrophone. 

The  measured  decibel  gain,  (47)— (49) ,  of  transfer  function 

(61)  is  displayed  in  figure*  2  for  the  range  of  frequencies  from 

30  to  500  Hertz,  on  a  logarithmic  frequency  abscissa.  Also 

superposed  are  the  decibel  gain  responses  of  filters  with  1  or  2 

or  3  poles  at  the  origin,  which  plot  as  straight  lines  on  this 

type  of  paper.  This  information  is  required  for  determining  the 

behavior  of  the  filter  from  30  Hertz  down  to  f  -  0  and  is 

necessary  because  the  Hilbert  and  Fourier  procedures  both  require 

knowledge  of  the  complete  attenuation  (or  gain)  for  all 

frequencies,  in  order  to  determine  the  value  of  the  corresponding 

minimum-phase  shift  at  just  one  frequency.  It  may  be  reasonably 

concluded  from  the  fits  in  figure  2  that  the  transducer  of 

* 

Figures  2  through  11  are  collected  at  the  end  of  this  section. 
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interest  here  has  a  double  zero  at  f  -  0. 

In  addition,  the  same  fitting  procedure  has  been  attempted  in 
the  neighborhood  of  500  Hertz  in  figure  2,  as  may  be  seen  by  the 
superposition  of  responses  for  filters  with  decays  corresponding 
to  0  or  1  or  2  or  3  poles  at  f  ■  «.  However,  the  situation  is 
rather  poor  at  this  upper  end  of  the  measured  frequency  range, 
because,  as  seen  in  figure  2,  the  transducer  has  not  yet 
developed  its  asymptotic  behavior  at  f  ■  500  Hertz.  This 
behavior  is  consistent  with  the  information  mentioned  above, 
which  describes  the  use  of  this  device  as  a  source  up  to  900 
Hertz.  Thus,  we  have  a  situation  where  we  have  insufficient 
measurements  to  fully  apply  the  theoretical  developments 
presented  earlier.  Nevertheless,  we  will  attempt  to  circumvent 
the  inadequacy  by  extrapolating  the  given  measurements  into  the 
frequency  range  above  500  Hertz  and  then  using  the  combination  of 
measured  and  extrapolated  gains  to  determine  the  minimum-phase 
response . 

PHILOSOPHY  OF  EXTRAPOLATION 

A  situation  of  frequent  occurrence  is  the  following.  We  have 
a  measured  residual  attenuation  ctjff),  but  it  is  available  only 
for  0  <  f^  <  f  <  fj;  see  (54)— (57).  We  presume  that  attenuation 
otjff)  is  even  about  f  -  0.  Call  this  total  frequency  range  of 
known  values,  K.  Denote  the  remainder  of  the  frequency  range, 
where  {^(f)  is  unknown,  by  U. 

We  want  to  evaluate  the  minimum-phase  corresponding  to  a-(f) , 
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namely 


02(f)  -  -  5(«2(£>)  -  -  H 


du 


<x2(u) 
f  -  u 


(62) 


Our  approach  is  to  extrapolate  a2(f)  beyond  K  into  the  unknown 
frequency  range  U.  Call  this  extrapolated  function  ®2e(f)»  it 
exists  for  all  f.  This  extrapolation  must  be  rather  close  to  the 
true  (unknown)  attenuation  a2 ( f )  in  U,  but  a2g(f)  need  not  agree 
with  a2  ( f )  inside  K.  In  particular,  a2e(f)  and  a2  ^  ^  should 
match  in  value  and  slope  at  the  boundaries  of  K. 

Then,  we  can  obtain  the  following  approximation  to  phase 
( 62 ) ,  namely 

a,(u)  ,  r  a2  (u) 

du 


e  (f)  e  .  I  [  du  _  I  f 

2a  n  J  f  -  u  r  J 


f  -  u 


U 


-H 


du 


«2(U)  - 


«2e(u) 


+  ® 


f  -  u 


-H 


du 


a2e(u) 

f  -  u 


K 


(63) 


The  first  (finite)  integral  in  (63)  is  done  numerically,  by 
employing  the  Fourier  procedure  presented  here.  The  second 
integral  in  (63)  is  actually  divergent  and  is  instead  replaced  by 
use  of  a  known  attenuation/minimum-phase  pair,  ^2e^^* 

The  key  to  this  procedure  is  a  shrewd  choice  for  the 
extrapolated  attenuation  a2e(f).  Several  candidates,  along  with 
the  corresponding  minimum-phase  functions,  are  listed  in  appendix 
D. 
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LAPLACE  TRANSFORM  NOTATION 

For  convenience  of  notation,  we  employ  here  the  Laplace 
transform  of  the  impulse  response,  namely 

+« 

L ( s )  -  J  dr  exp(-sr)  h(r)  ,  (64) 

0 

where  we  have  specifically  limited  consideration  to  causal 
filters.  The  connection  with  the  Fourier  transform  (1)  is 

H  (  f  )  «=  L(i2nf)  .  (65) 


EXAMPLE  A 


The  first  attempted  fit  to  the  measured  gain  in  figure  2  is 
by  means  of  filter 


L(s) 


c  s 


( s  +  a ) ( s  +  b) 


(66) 


with  constants  a  -  260,  b  -  330,  and  c  -  -  .55E8.  This  filter 
has  the  desired  double-order  zero  at  the  origin,  but  does  not 
decay  for  large  frequencies.  The  gain  of  (66)  is  superposed  on 
the  measured  gain  in  figure  3;  it  is  seen  that  the  constants  have 
been  chosen  to  give  a  fit  that  matches  in  value  and  slope  for 
small  frequencies  and  that  matches  the  measured  gain  value  at 
500  Hertz. 

The  difference  in  decibels  between  the  measured  gain  and  the 
fitted  gain  is  displayed  in  figure  4;  it  goes  to  zero  at  30  and 
500  Hertz  and  is  assumed  to  be  zero  outside  this  range.  This 
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assumption  is  not  likely  to  be  correct  for  f  greater  than  500 
Hertz,  but  it  is  necessary  in  order  to  proceed  with  the  numerical 
manipulations.  The  difference  in  attenuations,  e^ff)  of  (55),  is 
available  by  dividing  the  result  in  figure  4  by  -8.686;  see 
( 47 )-( 49  )  . 

The  residual  attenuation  {^(f)  is  subjected  to  the  cascaded 
Fourier  procedure  of  (46),  and  the  resultant  phase  (^(f)  is 
to  the  minimum-phase  8^ ( f )  corresponding  to  (66).  The  final 
total  phase  8(f)  is  shown  in  figure  5,  with  the  label  A&T, 
meaning  analytic  and  transform,  that  is,  81(f)  plus  82(f)* 
Superposed  on  this  figure  is  the  measured  phase,  with  the  label 
M&D,  meaning  measured  and  time-delay  adjusted.  Recall  in  the 
discussion  surrounding  (61)  that  there  is  an  unknown  time  delay, 
between  the  transducer  and  receiving  hydrophone,  included  in  the 
measurements  taken.  Accordingly,  a  selection  of  time  delay  was 
made  that  yielded  the  best  eyeball  fit  of  the  two  phases  over 
the  range  of  frequencies  from  0  to  400  Hertz  in  figure  5;  this 
corresponds  to  an  additive  linear  phase  function  of  frequency,  as 
indicated  by  example  (60).  The  time  delay  was  1.43  ms. 

The  agreement  between  the  minimum-phase  and  the  measured 
results  in  figure  5  allow  us  to  conclude  that  the  J15-1 
transducer  is  indeed  a  minimum-phase  filter,  at  least  over  the 
frequency  range  up  to  400  Hertz.  The  difference  between  the  two 
results  is  17°  at  500  Hertz,  which  is  significant.  However,  the 
reason  for  this  discrepancy  is  undoubtedly  due  to  the  fact  that 
(66)  is  not  the  correct  fit  for  f  >  500  Hertz,  because  (66)  has 
no  decay  for  large  frequencies. 
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EXAMPLE  B 

In  an  effort  to  find  a  better  phase  match,  another  fit  was 
also  tried,  namely  filter 

c  s2 

L(s)  -  - — - = - =-  ,  (67) 

(s  +  ao)[(s  +  a)z  +  bz] 

with  constants  ao  -  4000,  a  -  260,  b  -  400,  and  c  -  -  .275E12. 

The  measured  and  analytical  decibel  gains  are  plotted  in  figure 
6,  while  the  decibel  difference  is  plotted  in  figure  7.  The 
corresponding  two  phase  plots,  obtained  by  an  identical  procedure 
to  that  described  in  example  A  above,  are  presented  in  figure  8. 
Now,  the  difference  in  the  two  phase  curves  at  500  Hertz  has 
decreased,  but  only  slightly,  to  14°.  Apparently,  the  unmeasured 
decibel  gain,  in  the  frequency  range  above  500  Hertz,  is  causing 
inaccurate  calculations  of  the  minimum-phase  in  the  region  just 
below  500  Hertz,  due  to  our  inability  to  correctly  extrapolate, 
by  means  of  (66)  and  (67),  to  what  the  filter  gain  truly  was  in 
that  frequency  range.  This  supposition  is  consistent  with  the 
observation  that  the  minimum-phase  at  a  particular  frequency  is 
largely  governed  by  the  (rate  of  change  of  the)  attenuation  in 
the  neighborhood  of  that  frequency  (2;  page  345],  The  agreement 
in  phase  results  for  the  lower  frequencies  comes  about  because 
errors  in  gain  measurements  above  500  Hertz  have  a  much  reduced 
effect  on  the  calculated  phase  at  low  frequencies. 
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EXAMPLE  C 

In  an  attempt  to  justify  this  conjecture,  an  estimate  of  the 
unmeasured  gain  in  the  frequency  range  from  500  to  900  Hertz  was 
made  and  is  illustrated  in  figure  9.  A  droop  of  7  dB,  centered 
at  565  Hertz,  has  been  added  and  is  annotated  by  the  phrase 
"augmented".  The  fit  is  again  (66),  with  the  same  constants  as 
used  for  example  A,  and  is  superposed  in  the  figure. 

The  two  phase  curves  are  illustrated  in  figure  10.  Now,  the 
discrepancy  between  the  two  results  is  negligible  (within 
measurement  error)  all  the  way  up  to  500  Hertz,  the  maximum 
frequency  at  which  the  phase  was  measured.  Thus,  we  feel 
justified  in  concluding  that  the  device  under  investigation  is 
indeed  a  minimum-phase  filter,  at  least  over  the  measured 
frequency  range  up  to  500  Hertz. 

LIMITED  FREQUENCY  RANGE 

It  has  been  stated  above  that  the  measured  filter  appears  to 
be  minimum-phase  in  a  particular  frequency  range.  Strictly,  this 
is  not  a  valid  concept;  but  it  is  necessary  to  allow  for  it  in 
practice,  where  filter  responses  cannot  possibly  be  measured  for 
all  frequencies.  For  example,  suppose  that  the  transfer  function 
H(f)  has  a  collection  of  poles  and  zeros  in  the  upper-half 
f-plane,  all  fairly  near  the  origin  f  -  0.  In  addition,  let  H(f) 
have  a  pole-zero  pair  far  away  from  the  origin,  but  symmetrically 
located  about  the  real  f  axis,  so  as  not  to  affect  the  gain  or 
attenuation;  see  the  pair  near  f  ■  f~  in  figure  11. 
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Obviously,  the  filter  in  figure  11  cannot  be  minimum-phase, 
because  it  has  a  zero  in  the  lower-half  f-plane.  Yet,  its 
measured  phase,  for  frequencies  less  than  f^,  would  be 
indistinguishable  from  that  of  the  minimum-phase  filter  that 
does  not  contain  that  extra  pair.  Thus,  we  would  reasonably 
conclude,  upon  the  basis  of  the  measurements  made,  that  the 
filter  is  "minimum-phase  for  f  <  f ^ . "  Furthermore,  this  is  a 
practically  useful  concept  because  compensation  of  the  filter  in 
this  same  frequency  range  is  certainly  possible  and  allowable. 

In  other  words,  measurement  in  a  limited  frequency  range  only 
allows  us  to  make  conclusions  in  that  same  range;  in  fact,  the 
situation  is  slightly  worse  than  that,  because  the  edges  of  the 
range  may  also  be  open  to  question. 
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f  (Hertz) 

Figure  E.  Measured  Filter  Gain 


f  (Hertz) 


Figure  3.  Fitted  Gain  For  Example  fl 
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Figure  4.  Decibel  Difference  for  Example  fl 
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Figure  5.  Measured  and  Transformed  Phases  for  Example  fl 
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measured 
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Figure  G.  Fitted  Gain  for  Example  B 
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Figure  7.  Decibel  Difference  for  Example  B 
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igure  8.  Measured  and  Transformed  Phases  for  Example  B 
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Figure  9.  Fitted  Gain  for  Example  C 
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Figure  10.  Measured  and  Transformed  Phases  for  Example  C 
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SUMMARY 

For  a  minimum-phase  filter,  the  phase  shift  8(f)  can  be  found 
from  the  attenuation  a(f)  by  means  of  two  cascaded  fast  Fourier 
transforms,  once  the  logarithmic  singularities  in  o(f)  have  been 
subtracted  out  and  handled  analytically.  A  partial  accuracy 
check  is  automatically  built  into  the  procedure,  because  the  real 
part  of  the  output  should  agree  with  the  given  input;  the 
imaginary  part  of  the  output  is  the  desired  minimum-phase  result. 
This  Fourier  approach  yields  the  entire  phase  curve  for  all 
frequencies,  not  just  a  point-by-point  output,  as  a  Hilbert 
transform  numerical  integration  would  give. 

In  order  to  use  this  procedure,  the  attenuation  must  be 
measured  for  all  frequencies,  or  at  least  for  large  enough  and 
small  enough  frequencies  that  the  asymptotic  behavior  is  well 
developed  and  obvious.  A  plot  of  the  attenuation  (or  decibel 
gain)  on  a  logarithmic  frequency  abscissa  is  recommended  for  this 
purpose,  because  the  filter  magnitude  characteristic  should 
approach  a  straight  line  with  a  decay  equal  to  a  multiple  of 
6  dB/octave  in  the  neighborhood  of  zero  and  infinite  frequencies. 
Failure  to  make  a  complete  set  of  measurements  will  lead  to  the 
need  for  extrapolation  and  the  attendant  errors  that  can  occur 
with  such  a  procedure,  as  illustrated  here.  Furthermore, 
statements  about  the  minimum-phase  behavior  of  a  particular 
filter  can  only  be  made  (with)in  that  same  frequency  range. 
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APPENDIX  A.  PRINCIPAL  VALUE  INTEGRAL  EVALUATION 


Through  a  change  of  variable,  a  principal  value  integral  can 
be  put  in  the  form 

b 

I  ■  J  dt  ,  where  g  { 0 )  t  0  .  (A-l) 

-b 


Limit  b  can  be  finite  or  infinite.  (For  example,  (8)  fits  this 
form  when  we  let  g(t)  -  G(x-t)/n.)  Although  (A-l)  is  a  principal 
value  integral,  it  can  be  expressed  as  (ordinary  integrals) 


I 


9o(t) 

t 


b 

2J dt 


V1’ 

t 


b 

|  ^  [ g( t)  -  g( -t ) ] 


(A-2) 


-b  0  0 

where  gQ(t)  is  the  odd  part  of  g(t);  see  definition  (5).  This 
form  can  be  used  for  numerical  evaluation  whether  b  is  finite  or 
not.  If  b  is  infinite,  the  integrand  of  the  last  integral  in 
(A-2)  maintains  the  same  decay  with  t  as  original  integral  (A-l). 
This  is  not  true  of  the  sometimes  recommended  alternative  form 


i  .  J  at  am  -  9(0)  _  (A_ 

-b 

which  decays  very  slowly  with  t,  although  it  is  finite  at  the 
origin  t  -  0.  However,  another  alternative  that  advantageously 
uses  thi6  subtraction  device  is  given  later  in  (A-ll). 

A  simple  example  of  (A-l)-(A-2),  for  b  finite,  is  furnished 
by  the  integral 
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b 

dt  expU) 
-b 


sinh(t) 

at  — j —  . 


(A-4) 


the  latter  of  which  has  a  well-behaved  integrand  at  t  -  0. 


DERIVATIVE  EVALUATION 

In  general,  the  last  integrand  in  (A-2)  behaves  as 

g(t)  -  g(-t)  _  2  g,(0)  as  t  ■*  0  .  (A-5) 

Therefore,  in  order  to  use  (A-2),  it  is  necessary  to  have  g'(0). 
If  all  we  can  easily  evaluate  is  g(t),  and  not  its  derivative 
g'(0),  a  good  approximation  is  available  through  the  following 
device.  We  know  that  g'(0)  is  approximated  by 

g(e)  -  g(-e)  for  small  e  #  (a-6) 

Z  t 

However,  if  e  is  too  large,  this  is  a  poor  approximation,  whereas 
if  e  is  too  small,  round-off  errors  cause  numerical  stability 
problems.  But  we  know  that 

9(  ~  2c3(~C)  “  9'(°)  +  |  g' ' ' ( 0)  e2  +  0(e4)  as  e  0  .  (A-7) 

So,  letting  F(e)  be  the  left-hand  side  of  (A-7),  we  have,  to 
second  order, 

F(e)  -  Aq  +  Ax  e2  } 

*  \  where  A  and  A.  are  unknown  .  (A-8) 

F( e/2 )  -  Aq  +  Ax  eZ/4j 
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The  desired  unknown  follows  easily  from  (A-8)  as 

A0  -  4  r(t^>  -  M  .  ,.,01  . 


(A-9) 


This  procedure  is  an  extrapolation  to  the  limit;  it  uses  c/2  as 
the  smallest  argument  of  F. 

A  program  for  the  evaluation  of  g'(t)  at  general  t  is 
furnished  here  in  BASIC;  it  requires  specification  of  a  tolerance 
Tol  in  line  70  of  the  function  subroutine  FNDerivl. 


10 

INPUT  T 

20 

Deri- FNDerivl ( T) 

30 

PRINT  T, Deri  ! 

t,g' ( t) 

40 

END 

50 

! 

60 

DEF  FNDerivl(T)  ! 

~g'(t)  via  extrapolation 

70 

Tol-l.E-6  ! 

tolerance 

80 

E-.  2  ! 

epsilon  (start) 

90 

E-E* . 5 

100  V1-V2 

110  V2- ( FNG ( T+E ) -FNG ( T-E ) )/(2.*E) 

120  V-V2  + ( V2-V1 )/3 . 

130  IF  ABS( V2/V-1 . ) >Tol  THEN  90 
140  RETURN  V 
150  FNEND 
160  1 

170  DEF  FNG(T) 

180  RETURN  EXP(T)  !  example  exp(t) 
190  FNEND 


An  application  of  this  program  to  the  exp(t)  example  in  line  180, 
at  argument  t  -  1.1,  yielded  an  error  of  -7.8E-13. 

If  we  instead  kept  terms  to  fourth  order  in  (A-7),  an 
extension  to  (A-8)  yields  approximation 


9-(0)  «  ji[64  r(f)  -  20  r(§)  ♦  ru>]  . 

This  procedure  uses  e/4  as  the  smallest  argument  of  F. 


( A-10 ) 
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AN  ALTERNATIVE  SUBTRACTION  PROCEDURE 
We  now  express  (A-l)  in  the  form 

I  -  j  dt  2lLL  -  J  dt  2llL  +  J  dt  2lli  '  (A- 

-b  -a  R 

where  limit  a  is  chosen  for  convenience  and  R  is  the  union 
(-b,-a)  U  (a,b).  Then,  as  done  in  (A-3), 

I  -  J  dt  9(t)  ;  9(°)  +  j  dt  alii  .  (A. 

-a  R 

These  are  both  ordinary  integrals  now.  The  first  integrand  is 
finite  at  t  -  0,  with  value  g'(0),  while  the  second  integrand 
maintains  its  original  decay  as  x  -*  ±b. 


SECOND  DERIVATIVE  EVALUATION 

The  procedure  presented  in  (A-5)-(A-9),  for  the  approximate 
evaluation  of  first  derivative  g'(0),  can  be  extended  to  the 
second  derivative  g"(0)  as  follows.  We  know  that 


q(e)  ♦  q(-e) 


-  g(  0 ) 


g" ( 0 )  e2  +  0(e4) 


as  e  0 


(A-13) 


Therefore , 


g( e)  +  q(-e )  -  2q(0) 


g" ( 0  )  + 


0 (  c2 ) 


(A-l 4 ) 
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Letting  D(e)  be  the  left-hand  side  of  (A-14),  we  have,  to  second 
order , 


D(e) 

D(e/2) 


Bo  +  B1 


Bo  +  B1 


where  B  and  B,  are  unknown 

O  1 


(A-15) 


The  desired  solution  is 


4  D( c/2)  -  D( c) 


g"(0) 


(A- 16 ) 


This  is  an  extrapolation  to  the  limit;  it  uses  c/2  as  the 
smallest  argument  of  D.  A  program  for  the  evaluation  of  g"(t)  at 
general  t  is  given  below  in  BASIC;  it  requires  specification  of  a 
tolerance  Tol  in  line  70  of  the  function  subroutine  FNDeriv2. 


10 

INPUT  T 

20 

Der2-FNDeriv2(T) 

30 

PRINT  T , Der 2  ! 

t,g"(t) 

40 

END 

50 

j 

60 

DEF  FNDe r i v2 ( T )  ! 

~g"(t)  via  extrapolation 

70 

Tol-l.E-6  ! 

tolerance 

80 

E-.2  ! 

epsilon  (start) 

90 

G2-2 . *FNG( T ) 

100 

E-E* . 5 

110 

V1-V2 

120 

V2- ( FNG ( T+E ) +FNG ( T- 

E ) -G2 )/( E*E ) 

130 

V-V2+ ( V2-V1 )/3 . 

140 

IF  ABS ( V2/V-1 . ) >Tol 

THEN  100 

150 

RETURN  V 

160 

FNEND 

170 

! 

180 

DEF  FNG(T) 

190 

RETURN  EXP(T)  ! 

example  exp(t) 

200 

FNEND 

An  application  of  this  program  to  the  exp(t)  example  in  line  190, 
at  argument  1.1,  yielded  an  error  of  1.6E-11. 
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APPENDIX  B.  FOURIER  TRANSFORM  OF  GENERALIZED  FUNCTION 

We  are  interested  in  finding  the  Fourier  transform  of  the 
generalized  function 

exp(-ax)  u(t)  ^  a  >  0  ,  (B-l) 

where  U(x)  is  the  unit  step  function.  Letting  u  -  2nf,  the 
integral  of  interest  is 

I  -  J  ^  exp(-ax)  U(x)  exp(-iux)  - 
■  J  —  (exp(-ax)  -1+1]  U(x)  exp(-iux)  - 

+  CB 

■  -  J  ^  (1  -  exp(-ax)]  exp(-iux)  +  J  —  U(x)  exp(-iux)  - 
0 


-  ln(s-nrJfi)  -  [‘!  ‘9"  (2?)  +  ln 


0) 
2  ji 


+  C' 


( B-2  ) 


-  ln(a  +  iu)  +  ln(iu)  -  ij  sgn(u)  -  ln|u|  +  ln(2n)  -  C' .  (B-3) 


In  (B-2),  we  used  [4;  page  334,  3.434  2]  and  (6;  page  43,  row  3, 
column  3,  with  m  ■  1],  But  since 


ln ( iu) 


in/2  +  ln|u|  for  u  >  0 


l-iii/2  +  ln|u|  for 


u  >  O'! 
u  <  oj 


+  i2  Jin 


ij  sgn(u)  +  ln|u|  +  i2nn  ,  n  integer 


( B-4 ) 


45 


TR  8667 


we  can  express  ( B— 3 )  as 


I  -  -  ln(a  +  ico)  +  C  ,  where  C  -  ln(2n)  -  C'  +  i2nn  .  (B-5) 


Thus,  we  have  the  Fourier  transform  pair 


exp(  aT)  U(t)  « - -  ln(  a  +  i2nf)  +  C  , 


(B-6) 


where  C  is  an  arbitrary  constant.  The  reason  for  the  presence  of 
C  is  that  the  generalized  function  ^  U(t)  is  indeterminate  within 
an  additive  arbitrary  multiple  of  the  delta  function  5(t). 

For  the  example  in  (33)  of  H(f)  -  1/(1  +  i2nf),  we  have 
Q(f)  -  ln(l  +  i2nf).  Application  of  pair  (B-6),  with  a  ■  1,  to 
(39)  then  yields  causal  function 


q(  t) 


exp(-T) 


U(t) 


( B-7  ) 
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APPENDIX  C.  HILBERT  TRANSFORM  MANIPULATION 


It  was  noted  below  (45)  that  the  Hilbert  transform  of 
attenuation  a(f)  encounters  integrals  with  logarithmic  infinities 
and  must  be  handled  more  carefully.  This  problem  is  treated  in 
[3;  pages  206-8],  by  dividing  the  attenuation  by  a  factor  that  is 
quadratic  in  f,  rather  than  linear.  In  current  notation,  that 
result  is  [3;  (10-67)] 


(C-l) 


—  00 


If  we  utilize  the  property  employed  in  [3;  page  208,  line  2], 
namely  that  attenuation  a(f)  is  even,  we  can  develop  (C-l)  as 


+  CD 


1 


n 


du  a(u) 


( C-2  ) 


0 


( C—  3 ) 


0 


( C—  4  ) 
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The  step  leading  from  (C-2)  to  (C-3)  presumes  that  both  of 
the  latter  integrals  converge  separately,  which  need  not  be  the 
case  for  attenuations  a(f);  this  is  the  reason  for  the  quadratic 
denominator  adopted  in  (C-l),  which  guaranteed  convergence  of 
that  integral. 

Rather  than  using  Hilbert  transforms  and  having  to  employ  the 
method  of  (C-l),  we  have  resorted  instead  to  the  use  of  Fourier 
transforms,  as  outlined  in  (46).  Of  course,  a  similar  problem 
arises  there,  as  mentioned  in  the  sequel  to  (51).  The  method  of 
circumventing  the  difficulty,  in  the  Fourier  approach,  is  to 
subtract  out  the  singularities  and  handle  them  analytically,  as 
described  in  (54)— (57) . 

The  justification  of  this  procedure,  using  modified  Hilbert 
transform  (C-l)  as  a  starting  point,  is  as  follows.  Express 
given  attenuation  a(f)  in  two  parts,  as  in  (54),  where  residue 
<*2  ( f )  has  a  convergent  Hilbert  transform  integral 


for  all  f  . 


(C-5) 


The  phase  shift  8(f)  corresponding  to  attenuation  a(f)  is  then 
given  by  sum  (56),  where,  following  (C-l), 


ex(f ) 


( C-6 ) 


and  82(f)  is  available  as  the  negative  of  (C-5).  The  proof  of 
this  last  claim  follows  immediately  from  the  derivation  in 
(C-l)-(C-4)  if  we  replace  a(f)  and  8(f)  everywhere  by  o^f)  and 
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f$2  ( f )  •  respectively.  This  is  legitimate  because  the  existence  of 
(C-5)  for  residual  attenuation  ajff)  now  a^ovs  the  separation 
into  two  convergent  integrals,  as  done  in  (C-3). 

We  do  not  actually  use  (C-5)  or  (C-6).  Instead,  (C-6)  is 
accomplished  by  using  known  closed  form  attenuation/minimum-phase 
pairs  for  o^(f)  and  8^(f),  while  (C-5)  is  replaced  by  the  Fourier 
approach  given  in  (46),  with  0.2(f)  and  substituted  for  a(f) 

and  8(f),  respectively.  The  inverse  Fourier  transform  integral 
in  the  top  line  of  (46),  but  now  in  terms  of  aj(f )»  is 
conve  rgent . 

(For  interest,  an  example  of  the  application  of  (C-6)  is 
afforded  by  attenuation-phase  pair  (51).  This  fact  is 
immediately  verified  by  use  of  (4;  4.295  8].) 
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APPENDIX  D.  EXAMPLES  OF  ATTENUATION/MINIMUM-PHASE  PAIRS 

In  this  appendix,  we  list  a  few  attenuation/minimum-phase 
pairs  that  can  be  used  in  the  subtraction  procedure  presented  in 
( 54  )  —  (  57)  to  eliminate  the  divergent  integrands  encountered.  For 
convenience  of  notation,  we  employ  the  Laplace  transform  of  the 
impulse  response,  namely 

+® 

L(s)  -  J  dr  exp(-sx)  h(r)  ,  (D-l) 

0 

where  we  have  specifically  limited  consideration  to  causal 
filters.  The  connection  with  the  Fourier  transform  (1)  is 

H ( f )  -  L(i2nf)  .  ( D-2 ) 

In  the  following,  a,  b,  and  c  are  real  positive  constants,  and 
o)  -  2  nf  . 

EXAMPLE  1 : 


L(s) 


s  +  a  ' 


o(f )  -  j  ln(a^  +  oi2)  -  ln(c)  ,  3(f)  ■  arctan(  <o/a )  .  (D-3) 


In  the  limit  as  a  ■*  0+, 


o(f)  -  In | w |  -  ln(c)  ,  3(f)  -  j  sgn(w)  . 


(D-4) 
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EXAMPLE  2: 


L(S) 


c  s 


s  +  a 


o(f)  -  |  ln( a2  +  to2)  -  In  |  co|  -  ln(c)  , 


8(f)  -  arctan(  co/a )  -  j  sgn(co) 


EXAMPLE  3 


Ms) 


c  s 


(  s  +  a  j  (  s  +  b )  ' 


o(f)  -  |  ln(  a2  +  to2)  +  |  ln(  b2  +  o>2 )  -  In  |  co| 


8(f)  -  arctan(co/a)  +  arctan(co/b)  -  y  sgn(co) 

u 

This  attenuation  reaches  a  minimum  at  w  -  (ab)  ,  at 
the  phase  goes  through  zero. 

EXAMPLE  4: 


Ms) 


2  2  ' 
(s  +  a r  +  o 


o(f)  -  |  In  [( a2  +  (to  +  b)2]  +  |  In  [a2  +  (to  -  b)2] 


8(f)  -  arctan  ~  — j  +  arctan  • 


( D-5 ) 

ln( c )  , 

(D-6) 

which  point 

-  ln( c )  , 

(D-7) 
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APPENDIX  E.  NUMERICAL  EVALUATION  OF  (46) 

We  repeat  here  the  cascaded  Fourier  transform  operations 
listed  in  ( 46 ) : 

3(T)  «  F_1{«(f)}  ,  ( E— 1 ) 

q(x)  -  2  3(x)  U(x)  ,  (E-2) 

<%(f)  +  i  8(f)  -  F{q(x)}  .  ( E-3 ) 

We  limit  consideration  to  the  case  where  attenuation  a(f)  is 
even,  which  is  the  typical  practical  situation.  Also,  we  weight 
the  inverse  Fourier  transform  in  ( E— 1 )  by  real  symmetric  window 
W(f),  which  is  zero  for  |f|  >  MA.  We  then  get  approximation 

+ «° 

<jQ(x)  *  J  df  exp(i2nfx)  o(f)  W(f)  - 
•  0 

+  oo 

-  2  Re  J  df  exp(-i2nfx)  o(f)  W(f)  - 

0 

MA 

-  2  Re  J  df  exp(-i2nfx)  a(f)  W(f)  ■ 

0 

M 

b  2  Re  )  |  s  A  exp( -i2imAx)  a(nA)  W(nA)  ■  ,  ( E-4 ) 

n-0  n 

where  we  sample  in  frequency  f  with  increment  A.  We  also  use 

some  integration  rule  like  trapezoidal  or  Simpson;  for  example, 

the  trapezoidal  rule  has  s  ■  1,  except  for  s  -  s„  -  1/2. 

n  on 
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The  approximation  (jb(x),  defined  by  the  bottom  line  of  (E-4), 
has  period  1 /t  in  t.  Therefore,  we  compute  it  at  the  points 

T  -  jtx  for  0  <  m  <  N  -  1  ,  (E-5) 

N  A 

which  cover  a  full  period  of  <jb(x).  There  follows 
(  \  M 

ShImtI  “  2A  Re  )  |  s  exp(-i2nnra/N)  a(nb)  W(n&)  ,  (E-6) 

dvng;  n_0  n 

which  is  an  N-size  fast  Fourier  transform  of  M  +  1  data  points. 
Any  surplus  points  can  be  collapsed,  if  desired,  without  loss  of 
accuracy;  see  [7;  pages  4-5],  for  example. 

Operations  ( E-2 )  and  (E-3)  can  be  combined  to  read 

+® 

Q(f)  -  a(f)  +  i  8(f)  -  2  J  dx  exp(-i2nfx)  <j(x)  .  (E-7) 

0 

Because  all  we  have  available  is  approximation  gb(x)  from  (E-4), 
we  adopt  the  following  approximation  to  Q(f),  based  on  (E-7): 

+<*> 

Qa(f)  ■  2  J  dx  exp(-i2nfx)  <jb(x)  - 
0 

.S/t 

=  2  J  dx  exp(-i2nfx)  <jb(x)  -  (E-8) 

0 


2  "m  S5  e*p(-i2"£N?)  SbfcS)  *  °b(£; 

m-u 


(E-9) 
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where  wm  is  an  integration  weight.  The  integral  in  (E-8)  was 
limited  to  .5/A  in  t,  because  approximation  ^(t)  in  (E-4)  is 
only  available  up  to  that  limit  without  aliasing. 

The  period  of  the  final  approximation  Qjj(f)  in  (E-9)  is  NA  in 
f.  Therefore,  we  limit  its  computation  to  the  values 


for  0  <  n  <  N-l  .  ( E-10 ) 


This  can  be  accomplished  as  an  N-size  fast  Fourier  transform  of 
N/2  +  1  data  points.  The  final  approximation  to  desired  phase 
(3(f)  in  ( E-7 )  is  available  as  the  imaginary  part  of  (E-10),  at 
frequencies  f  -  nA.  In  addition,  the  real  part  of  (E-10)  should 
be  in  very  good  agreement  with  specified  attenuation  values 
(o(nA)  W(nA)}  used  in  (E-6);  this  serves  as  an  accuracy  check  on 
the  complete  procedure.  Equations  (E-6)  and  (E-10)  are  the  final 
results.  Strictly,  (E-6)  should  be  applied  only  to  the  residual 
attenuation  «2 ( f )  defined  in  (55);  then  (E-10)  furnishes  an 
approximation  to  t^tf)  +  i  PjCf)*  A  program  in  BASIC  for  the 
Hewlett  Packard  9000  computer,  for  the  procedure  given  above,  is 
presented  below. 
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10  ! 

NIJSC  TR  8667,  FOURIER  PROCEDURE  APPLIED 

20  ! 

TO  REAL  EVEN  FUNCTION  OF  FREQUENCY 

30 

De 1 1  af  =  5 .  ! 

SAMPLING  INCREMENT  IN  FREQUENCY 

40 

Fmax=900.  ! 

MAXIMUM  FREQUENCY 

50 

N= 16384  ! 

SIZE  OF  FFT 

69 

fl=260 .  ! 

FILTER  PARAMETERS 

70 

Btt330 .  ! 

FOR 

80 

C  =  - . 55E8  ! 

EXAMPLE  C 

90 

COM  fl,B,C 

100 

REDIM  Cos<0s  N/4) , X<0: N-l ) , Y(0:N-1 > 

110 

DIM  Cos< 4096) , X< 16384), Y< 16384), Real 

euen(25000),Phase(6s 100) 

120 

DOUBLE  N,M,Ns,Ms,N2,M2  ! 

INTEGERS 

130 

T=2. *PI/N 

140 

FOR  Ns*0  TO  N/4 

150 

Cos(Ns)=COS(T*Ns)  ! 

QUARTER-COSINE  TABLE 

160 

NEXT  Ns 

170 

M=Fmax/De 1 t  af 

180 

REDIM  Realeuen<0:M) 

190 

CALL  Input  real  even < De 1 t af , F max , Re al even <*) )  !  RESIDUAL 

200 

MAT  X= ( 0 . )  ! 

ATTENUATION  ALPHA2 

210 

MAT  Y  =  (  0 . ) 

220 

XC0>=.5*Realeuen<0> 

230 

Ms  =  M  MODULO  N 

240 

X(Ms)=. 5*Real euen(M) 

250 

FOR  Ns- 1  TO  M- 1 

260 

Ms=Ns  MODULO  N  ! 

COLLAPSING 

270 

X(Ms)=X(Ms)+Realeuen(Ns) 

280 

NEXT  Ns 

290 

CALL  Fftl4(N,Cos(*),X(*),Y<*))  ! 

FOURIER  TRANSFORM 

300 

N  2  =  N  /  2  ! 

INTO  TIME  DOMAIN 

310 

G  I N I T 

320 

PLOTTER  IS  "GRAPHICS" 

330 

GRAPHICS  ON 

340 

WINDOW  -N2,N2,-6,2 

350 

LINE  TYPE  3 

360 

GRID  N/8, 1 

370 

PRINT  "FOURIER  TRANSFORM  (TIME  DOMAIN) " 

380 

FOR  Ns=-N2  TO  N2 

390 

Ms*Ns  MODULO  N 

400 

PLOT  Ns,LGT(ABS(X<Ms))+l.E-99)  ! 

TIME  DOMAIN  FUNCTION 

410 

NEXT  Ns 

420 

PENUP 

430 

PAUSE 

440 

MAT  Y* <0 • ) 

450 

T=4./N  ! 

2  De 1 1 af  *  2  ✓  <N  Deltaf) 

460 

FOR  Ms*0  TO  N2 

470 

X  <  Ms ) =X<  Ms  >  *T  ! 

DOUBLE  FOR  POSITIVE  TIME 

480 

NEXT  Ms 

490 

X<0>«X<0>*. 5 

500 

X < N2 >  sX< N2  > * • 5 

510 

FOR  Ms*N2+ 1  TO  N-l 

520 

X  <  Ms ) =0 .  ! 

ZERO  FOR  NEGATIVE  TIME 

530 

NEXT  Ms 

540 

CALL  Fft  14<N,Cos<0,X<*>,  Y<*>>  ! 

FOURIER  TRANSFORM 

550 

M2=M*2  ! 

INTO  FREQUENCY  DOMAIN 

56 


560 

570 

580 

590 

600 

610 

620 

630 

640 

650 

660 

670 

680 

690 

700 

710 

720 

730 

740 

750 

760 

770 

780 

790 

800 

810 

820 

830 

840 

850 

860 

870 

880 

890 

900 

910 

920 

930 

940 

950 

960 

970 

980 

990 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 
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GCLEflR 

WINDOW  0 , M2 , - 1 , 1 
LINE  TYPE  3 
GRID  N/16, . 2 

PRINT  "ORIGINAL  INPUT  (FREQUENCY  DOMAIN)" 

FOR  Ns*0  TO  M  I N  (  M  , N 2 ) 

PLOT  Ns , Real even( Ns  >  !  ORIGINAL  INPUT 

NEXT  Ns 

PENUP 

PAUSE 

LINE  TYPE  1 
FOR  Ns*0  TO  M2 

PLOT  Ns , X (Ns )  !  F-T-F  APPROX  I  MAT  I  ON 

NEXT  Ns 

PENUP 

PAUSE 

DATA  -38.6,-48.2,-54.8,-60.4,-76.2,-82. 1,-94.5,-103.8,-109.1,-117.1 

DATA  -124. 1,-134.0,-143. 1,-152.9,-163. 1,-172.4, 179. 1, 171. 1, 164.2, 157.9 

DATA  152.8,147.1,142.8,135.8,131.9,128.7,122.8,118.7,115.1,110.6 

DATA  105.9,103.4,102.8,99.9,98.6,93.8,93.1,91.2,89.6,89.5 

DATA  89.6,89.6,89.2,88. 1,85.6,84.5,82.0,81. 1,79.0,74.7 

DATA  71.4,66.5,61.3,55.1,48.1,41.6,34.0,29.3,22.0,16.1 

DATA  12. 2, 5. 7, 2. 4, -3. 1,-6. 5, -11. 3, -16. 2, -2 1 . 2 , -25 . 7 , -29 . 7 

DATA  -33. 4, -37. 0, -40. 7, -43. 5, -47. 0, -49. 5, -51 . 6, -54. 1 , -56.2,-59.4 

DATA  -61 . 0, -62. 4, -64. 2, -66. 7, -68. 7, -71 . 4, -74. 6, -78. 1 , -8 1 . 4 , -83 . 8 

DATA  -88.7,-91.3,-95.0,-98.7,-103. 1 

READ  Phase ( * )  !  MEASURED  PHASE  IN  DEGREES 

FOR  Ns=22  TO  100 

Phase(Ns)=Phase(Ns)-360.  !  UN-WRAPPING  OF  PHASE 

NEXT  Ns 

MAT  Phase  =  Phase* ( -P I / 1 80  . )  !  MEASURED  PHASE  IN  RADIANS 

T*=2.  *PI*Del  taf 
FOR  Nse0  TO  N2 
W  =  T  *Ns 


Phaseapp=ATN( (W-B)/A)+ATN( (W+B)/A> 
X(Ns )«Phaseapp+Y(Ns  >  ! 

NEXT  Ns  ! 

GCLEAR 

WINDOW  0,  1 80 , 0 , P I  *  1 . 25 
LINE  TYPE  1 
GRID  20, PI*. 25 

PRINT  "PHASE  (FREQUENCY  DOMAIN)" 
FOR  Ns“0  TO  180 

PLOT  Ns , X ( Ns )  ! 

NEXT  Ns 

PENUP 

LINE  TYPE  3 
FOR  Ns *6  TO  100 

PLOT  Ns, Phase(Ns)-Ns*.0448  ! 

NEXT  Ns  ! 

PENUP  ! 

PAUSE 
END 
! 


!  PHASE  BETA  1  OF  APPROX. 
CALCULATED  PHASE  IN  RADIANS: 
BETA  *  BET  A 1  +  BETA2 


PHASE  VIA  FOURIER  PROCEDURE 


MEASURED  PHASE  WITH 
TIME  DELAY  CORRECTION 
OF  1.43  MILLISECONDS 
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1100 

1110 

1120 

1130 

1140 

1150 

1160 

1  170 

1  180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 

1450 

1460 

1470 

1480 

1490 

1500 

1510 

1520 

1530 

1540 

1550 

1560 

1570 

1580 
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SUB  Fft 1 4  <  DOUBLE  N,  REAL  Co£.<*),X<*>  Y(*n  .  ^  o 
DOUBLE  Log2n,Nl,N2,N3,N4  J  K  i  Ee  ,H<‘2M4"163M»  0  SUBS 
DOUBLE  K  ^ie  9  '  f,  =  2,147,483,648 

IF  N=1  THEN  SUBEXIT  1  1 1 1 ,  1 1 2 , 1 1 3 , I  1 4 , L < 0 : 1 3 > 

IF  N>2  THEN  1220 
fi=X<0)+X< 1 > 

X<1 >=X<0)-X< 1 > 
x<0)«n 
fl=Y<0)+Y< 1 > 

Y  < 1 >  *Y <  0  > - Y< 1 > 

Y<0)=fl 

SUBEXIT 

fl  =  L0G<N)/L0G<2.  > 

Log2n=fl 

IF  RBS<fl-Log2n>< 1 . E-8  THEN  1270 

PRINT  "N  -".N.-1S  NOT  A  POWER  OF  2i  DISALLOWED. " 

N1  *N/"4 


N2=N 1+1 
N3«N2+1 
N4=N3+N 1 

FOR  11=1  TO  Log2n 
I2=2-'<Log2n-I 1  > 

13=2*12 
1 4  =  N  / 1 3 

FOR  15=1  TO  12 
I6=< 15-1 >*I4+1 
IF  I  6 <  =  N 2  THEN  1410 
fil=-Cos<N4-I6-l > 

R2=-Cos< I6-N1-1) 

GOTO  1430 
R1 =Cos< 16-1 > 
R2=-Cos<N3-I6-l > 

FOR  17=0  TO  N-I3  STEP  13 
18=17+15-1 
19=18+12 
T 1 *X  < 1 8 ) 

T2=X( 19) 

T3=Y< 18) 

T4=Y< 19) 

R3=T1-T2 
R4=T3-T4 
X <  I8)*T1  +  T2 
Y  < 1 8 ) =T3  +  T4 
X< I9)=fll*fl3-fi2*fl4 
Y< I9)=fll*fl4+fi2*fi3 
NEXT  17 
NEXT  15 
NEXT  II 
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1590 

I 1 =Log2n+ 1 

1600 

FOR  12=1  TO  14 

1610 

LCI2-1 )=1 

1620 

IF  1 2  >Log2n  THEN  1640 

1630 

L  <  12—1  )=2/s(I  1-12) 

1640 

NEXT  12 

1650 

K  =  0 

1660 

FOR  11  =  1  TO  L < 1 3 ) 

1670 

FOR  12=11  TO  L < 1 2 )  STEP  L<13> 

1680 

FOR  13=12  TO  LC11)  STEP  L<12> 

1690 

FOR  14=13  TO  L ( 1 0 )  STEP  L <  1  1  > 

1700 

FOR  15=14  TO  L < 9 )  STEP  Ld0> 

1710 

FOR  16=15  TO  L C 8 )  STEP  L<9> 

1720 

FOR  17=16  TO  L C 7 )  STEP  L ( 8 ) 

1730 

FOR  18=17  TO  LC 6 )  STEP  L<7) 

1740 

FOR  19=18  TO  L  < 5 )  STEP  L<6> 

1750 

FOR  110=19  TO  L  C 4 )  STEP  L<5> 

1760 

FOR  111*110  TO  L  <  3 )  STEP  L<4> 

1770 

FOR  112=111  TO  LC2)  STEP  L<3> 

1780 

FOR  113=112  TO  Ld>  STEP  L<2> 

1790 

FOP  114=113  TO  L (0 )  STEP  Ld> 

1800 

J=I 14-1 

1810 

IF  K  > J  THEN  1880 

1820 

fl  =  X  oo 

1830 

X(K)=X< J) 

1840 

X< J)*R 

1850 

R=YCK) 

I860 

YCK)=Y< J) 

1870 

Y ( J ) =R 

1880 

K*K4  1 

1890 

NEXT  114 

1900 

NEXT  113 

1910 

NEXT  112 

1920 

NEXT  Ill 

1930 

NEXT  110 

1940 

NEXT  19 

1950 

NEXT  18 

I960 

NEXT  17 

1970 

NEXT  16 

1980 

NEXT  15 

1990 

NEXT  14 

2000 

NEXT  13 

2010 

NEXT  12 

2020 

NEXT  11 

2030 

SUBEND 

2040 

1 

59 


2050 

2060 

2070 

2080 

2090 

2100 

2110 

2120 

2130 

2140 

2150 

2160 

2170 

2180 

2190 

2200 

2210 

2220 

2230 

2240 

2250 

2260 

2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 
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SUB  Input_rea1_even<De1taf,Fmax,Rea1even<*)) 

DOUBLE  N s 

ALLOCATE  Db<6:180)  !  30:900  H2 

DATA  41.3,44.3,46.1,47.6,49.9,51.4,52.9,54.4,54.8,56.3 
DATA  57.0,57.4,57.9,58.6,59.0,59. 1,59.0,58.9,58.9,58. 8 
DATA  58.6,58.1,58.2,58.1,58.0,57.9,57.8,57.2,56.9,56.7 
DATA  56.6,56.4,56.3,56.2,55.7,55.6,55.4,55.0,54.9,55.2 
DATA  55. 2, 55. 7, 55. 7, 56. 1 , 56. 1 ,56. 6, 56. 9, 57. 5, 58. 3, 58. 6 
DATA  59.0,59.7,60.3,60.7,60.9,61.1,61.1,61.2,61.0,60.9 
DATA  60.7,60.6,60.4,60.2,60.0,59.9,59.6,59.4,59.3,58.7 
DATA  58.5,58.3,57.8,57.5,57.3,57.0,56.7,56.3,56.1,55.9 
DATA  55.7,55.5,55.6,55.6,55.4,55.3,55.4,55.3,55.6,55.3 
DATA  55.4,55.0,55.0,55.0,54.8 
REDIM  Db  <  6 :  100) 

READ  Db<*> 

MAT  Db=Db+<100.)  !  MEASURED  DB  GAIN 

REDIM  Db  <  6 : 180) 

FOR  Ns= 1 0 1  TO  180  !  AUGMENTED  DB  GAIN 

F  =  De 1 1  af  *Ns 
T  1  =  <  F-550 . )*. 04 
T2=  <  F-580 .  )*.  04 

Db<Ns)=154.8-5.*EXP<-Tl*Tl )  -  5 . *EXP<-T2*T2) 

NEXT  Ns 

MAT  Real evcna(0.  ) 

COM  A, B , C 
A2=A*A 
B2=B*B 
C2=C*C 

D1=<A2+B2)*<A2+B2) 

D2  =  2 . *  <  A2-B2 ) 

T*2 . *P I *De  1 1  af 
FOR  Ns  =  6  TO  180 


W  =  T  *  N  s 
W2CW*W 
W4=W2*W2 

PcC2*W4/<Dl+D2*U2+W4) 

Altenapp--.5*L0G<P) 

Atten=Db<Ns)/<-8.686) 

Real even<Ns)  =  At  t en-At t enapp 

NEXT  Ns 

SUBEND 


!  APPROX.  ATTEN.  ALPHA  1 
!  ATTENUATION  ALPHA 
!  RESIDUAL  ATTEN.  ALPHA2 
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